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Abstract 

The treatment of the time-independent Schrodinger equation in real space is an indispensable 
part of introductory quantum mechanics. In contrast, the Schrodinger equation in momentum 
space is an integral equation that is not readily amenable to an analytical solution, and is rarely 
taught. We present a numerical approach to the Schrodinger equation in momentum space. After 
a suitable discretization process, we obtain the Hamiltonian matrix and diagonalize it numerically. 
By considering a few examples, we show that this approach is ideal for exploring bound states in 
a localized potential, and complements the traditional (analytical or numerical) treatment of the 
Schrodinger equation in real space. 
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I. INTRODUCTION 



The treatment of the time- independent Schrodinger equation for a non-relativistic par- 
ticle of mass m is a prime element of quantum mechanics courses.-i^ The treatment of this 
second-order differential equation introduces students to the effect of boundary conditions 
on quantization^ and to the Sturm-Liouville problem.-i^ For only a few potentials V{x) can 
the Schrodinger equation 

be solved analytically, and the eigenvalue spectrum and the complete set of orthonormal 
eigenfunctions ipai^) be obtained. Introductory texts typically include a quantum well 
or step, a harmonic oscillator, a delta function potential, and various combinations.- The 
number of analytically solvable potentials is even smaller in higher dimensions. Notable 
exceptions are those with a central potential where rotational invariance allows us to obtain 
a second-order differential equation for the radial wavefunction in an effective potential that 
takes into account the centripetal barrier.-*^ Although such an equation is not, in general, 
analytically solvable, it is a significant improvement over the second-order partial differential 
equation. 

To explore the bound states in an arbitrary potential V{x), we can use the Wentzel- 
Kramers-Brillouin (WKB) approximation^ which provides a semiclassical picture of quan- 
tized energy eigenvalues. Another approach is to discretize Eq. ([1]) and obtain the matrix 
equation 



j=-N 



-I j=-N 

where Xj = jAx, Ax is the spacing between adjacent points along the discretized x-axis, and 
NAx = Xc denotes the spatial cutoff chosen such that X^.^ a where a is the characteristic 
length scale of the potential V{x). The tridiagonal second-derivative matrix has entries 
Dij = [Sij-i — 2Sij + Sij+i] /{Axy. In principle, as — oo and Ax — >• 0, the eigenvalues 
and eigenvectors of the (2A^ + 1) x (2A^ + 1) matrix Hij approach the spectrum of the 
original continuum problem. However, due to the diverging prefactor (Ax)"^ and the error 
in the matrix Dij at the end-points ±Xc, the matrix diagonalization approach does not lead 
to stable continuum results. Instead, the Numerov method has to be used to numerically 



obtain the eigenvalues and eigenfunctions.-i^ Another approach is to use the eigenfunction 
expansion method,- which results in a matrix equation for the expansion coefficients.- 

In this paper we show that the stability and convergence issues^ are circumvented by 
the Schrodinger equation in momentum space. In Sec. II we review the equation and the 
corresponding discretized-matrix eigenvalue problem. This method is ideal for localized 
potentials with a finite Fourier transform. We present the bound-state spectra for a few well- 
known potentials and compare them with analytical results whenever possible. In Sec. Illll we 
discuss the generalization of our approach to the Schrodinger equation in higher dimensions.- 
We conclude in Sec. HV] with a discussion and suggested problems. The method presented 
here is complementary to the standard differential-equation approach, can be explored in 
introductory quantum mechanics courses, and is accessible to junior or senior undergraduate 
students familiar with Matlab, Maple, Mathematica, or LAPACK. 



II. SCHRODINGER EQUATION IN MOMENTUM SPACE 

We start with the Fourier transform of the one- dimensional Schodinger equation, Eq. (fl]), 

f°° dp' 

ep^a(p) + J ~^')^"*^^') = (3) 

Here ipa{p) is the momentum-space wavefunction, which represents the probability amplitude 
for the particle to have momentum p, = /2m is the (non-relativistic) kinetic energy of 
the particle, and V{q) is the Fourier transform of the external potential V{x). We use the 
same symbol for the eigenvector \il)a) and the potential energy operator V in both real and 
momentum space; the fact that ipa{x) and ipaip) are different functions is understood.->i^ 
The integral equation ([3]) has been used to study the scattering problem^^!^ and the bound 
state in a 5-function potential.— 

To convert Eq. ([3]) into a form suitable for numerical exploration, we use Oq to denote 
the length-scale, define the momentum scale by h/ao, and use Eq = ff/2mal as the unit 
of energy. Because is arbitrary in the continuum limit, the spectrum with a potential 
characterized by depth Vq and range a, V{x) = Vof{x/a), is determined by the dimensionless 
parameter Vo/Ea = Vo/(^^/2ma^). For the numerical calculations it will be useful to choose 
Vo/Eq = 1 and vary a/ao to access various values of Vo/Ea- In terms of the dimensionless 
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variables Eq. ([3]) leads to a matrix equation after discretization: 



tpa{Un) = Hmn^aM = ^aVal^^m), (4) 



where a sum over the repeated index n is understood. Here Un = nAu = nApao/h = pnO-o/h 
is the dimensionless momentum, Am is the spacing between adjacent points along the u- 
axis, Ea = Ea/ Eq is the dimensionless eigenvalue, and V^n = V{Pm — Pn)^u/ {27rEQao) is 
the dimensionless potential along with the discrete measure Au/{2tt). Thus, the integral 
Schrodinger equation ([3]) has been recast as a matrix eigenvalue problem, Eq. (jlj). The 
eigenvalues and eigenvectors of the dimensionless Hamiltonian matrix ifmn are obtained 
using standard software packages. The results presented here were obtained by using Matlab 
and were verified by using LAPACK. This discretization process involves some computational 
subtleties that we discuss in the following, but provides an excellent way to study bound 
states in potentials that are localized in real space. 

The size of the Hamiltonian matrix Hmn is determined by the dimensionless ultraviolet 
momentum cutoff Uc = Pcdo/h, where is upper limit of the integration range in Eq. ([3]), 
and the dimensionless spacing Am. Note that 2TTao/Au and 2'n'ao/Uc impose the upper and 
lower limit respectively on the real-space size of the bound-state wavefunction. Thus, they 
need to be chosen for a given potential so as to obtain results that are valid in the continuum 
limit, Uc ^ oo and Au —>■ 0. If the external potential V{x) is even, the Hamiltonian Hmn is 
real. Therefore, the eigenfunctions ipa{p) have a definite parity and are real.- In this case it 
is sufficient to restrict ourselves to positive momenta m, n > 0. 

To demonstrate these considerations, we start with a quantum well with depth Vq and 
width a centered at the origin, V{x) = —Vo6{a — 2\x\) = V{—x) where 6{x) is the Heavy-side 
function.^- In this case the real Hamiltonian matrix is given by 



TT _ 2c _ Am 2aVo sin - Um)a/ao] _ ^ 
2n aoEo [(m„ - M„)a/aoJ 
The factor of 2 in the potential matrix elements arises from the restriction m,n > 0. The 

one-dimensional 5-function potential, V{x) = —Xi6{x), is obtained as a limiting case when 

a — > 0, Vq — i> oo with aVo = Ai. In this limit for an attractive potential, Ai > 0, the system 

has one exponentially bound state ipb{x) = ^/K,ex^p{—K\x\) with size = h^/mXi and 

energyi''^ E^ = —mX\/2h^. 

We obtain the eigenvalues and eigenvectors of the matrix Hmn for < \\i/EQaQ\ < 1 

with two different values of An ={0.01,0.005} and two different cutoffs Uc ={10,20}. The 



corresponding matrix dimension = Uc/ in these four cases varies from 1000 x 1000 to 
4000 X 4000. We find that the energy spectrum has one negative eigenvalue and the positive 
eigenvalues form a quadratic band representative of a free particle. (The positive-energy 
unbounded states are not accessible via the Numerov method.^) Figure [U^a) shows that 
the bound-state energy -Efe(Ai) matches the analytical result. Figure [T](b) shows a typical 
momentum-space wavefunction for the bound-state ipb{p) and a state with positive energy. 
As expected, we see that the bound-state wavefunction, ipbip) = 2(^k)^/^/(p^ + h^K^), is 
broad, whereas the positive-energy wavefunction is sharply peaked near a single momentum 
value. We check that the bound-state results are independent of Uc and Am. (The smallest 
momentum cutoff is chosen such that the contribution from momenta p > = UJi/aQ to the 
bound-state wavefunction is negligible.) Increasing Uc affects the eigenvalues and eigenvec- 
tors near the highest energy Ec/Eq = U^, whereas reducing Am sharpens the momentum- 
space eigenfunctions at positive energies. Thus, this numerical approach is particularly 
suited to study bound states, and may not handle the unbounded positive-energy states 
equally well. 

Next, we consider the problem of a deep quantum well Vo/Ea ^ 1. We diagonalize the 
matrix H„^n, Eq. ([5]), with Vq/Eq = 1, f/, = 300, Au = 0.1, and a/ao = {20,25}. Figure [2] 
shows that the numerically obtained spectrum of the bound-state energies measured from 
the bottom of the quantum well is quadratic. En = n^n'^S. This dependence is expected 
because for an infinite quantum well of size a the eigenvalues are given by En = n^Ti'^Ea. 
The prefactor £{Uc, A.u,Vq) ^ Ea, and a systematic exploration with increasing Uc and 
a, and decreasing Au shows that this discrepancy is due only to the discretization. We 
next consider an attractive Gaussian potential Vi(x) = — Vq exp(— x^/2a^). In this case a 
closed-form solution for the eigenvalues and eigenfunctions is unknown. The dimensionless 
Hamiltonian becomes 

H^^) _ ^ ^ exD [- -H(') (Q) 

Figure [3] shows the a-dependence of the magnitude of the ground-state energy E^ < Q 
obtained by using Vq/Eq = 1, Am = 0.01, and Uc = 30. The inset shows the ground state 
momentum-space wavefunction ipaip) for a/ao = {0.1,0.5}. As a/ao increases, the effective 
value of Vo/Ea increases. Thus, the ground state becomes more localized in real-space, and 
the spread of the wavefunction in momentum-space increases. 
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We emphasize that the bound-state eigenvalues and eigenf unctions obtained from the 
diagonahzation of the discrete matrix Hmn should be essentially independent of the cutoff 
Uc ^ 1 and the spacing Au <^ 1, to verify that they are valid in the continuum limit 
Uc oo, Au —>■ 0. Note that even in the limit Am ^ 0, a finite momentum-cutoff leads 
to a real-space potential that is not the same as the original one: 

^^^"^^ = /p ^ 
Thus the discretization parameters need to be so chosen that the difference between Vc{x) 
and V{x) is negligible. Two typical indicators that the continuum limit has not been reached 
are that some eigenvalues are lower than the depth of the potential well, Ea < —Vq, and the 
ground-state momentum-space wavefunction is linear, instead of quadratic, near p = 0. In 
one dimension we can choose the ground state wavefunction ipci^) to be positive, ^'^ so that, 
for an even potential the momentum-space wavefunction is parabolic at the origin, ipcip) — 
tpcip = 0) oc — p^. Therefore, a linearly varying ipcip) is a clear indication that the bound- 
state eigenvalues and eigenfunctions do not represent continuum results. In the following 
we show that the verification of the continuum limit is more subtle in two dimensions and 
requires a careful treatment. 



III. NUMERICAL APPROACH IN HIGHER DIMENSIONS 

For a particle in two or more dimensions, a naive discretization of the integral Schrodinger 
equation in Cartesian co-ordinate implies that the Hamiltonian matrix has a size ~ A^^ x A^^ 
where N = Uc/ Au is the number of discrete points along a single axis and D is the dimension. 
Thus, even in two dimensions, the parameters used in Sec. [Ill result in 10^ x 10^ or larger 
matrices that are impossible to treat numerically. For a central potential in two dimensions, 
the rotational invariance of the Hamiltonian implies that the angular momentum is a good 
quantum number and the eigenfunctions can be labeled by an integer angular momentum 
label i, ipa{p) = i^ae{p) exp{ii9p), where p = (p, 9p) is the two-dimensional momentum.-i^ii^ 
The Schrodinger equation for a given value of i becomes^ 

+ / ^^^^y(P^P" ^ - Op')e-'"^'^~'^'^^ai{p') = Ec^e^^eip) (8) 

where V{p,p'; 9p — 9pi) = V{\p — p'\) is the momentum-space potential and depends only on 
the angle between p and p' due to the central nature of the potential. The corresponding 
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dimensionless Hamiltonian matrix becomes Hmni^) = u^Smn + Vmn{(^), where the angular- 
averaged potential matrix is given by 

Here, < m„ < Uc denotes the magnitude of the dimensionless momentum and the matrix 
Hmn{^) has size N x N. Due to the prefactor Un^u from the two-dimensional area- 
element in polar co-ordinates, the Hamiltonian obeys Hnm{^) = {um/un)H^^{£). Thus, the 
discretized Hamiltonian matrix is not Hermitian with respect to transpose of the matrix 
plus complex conjugation. It is Hermitian with respect to the inner product defined via the 
two-dimensional measure. We will focus on £ = case because for a time-reversal invariant 
Hamiltonian, the ground state has zero angular momentum.-ii^ 

As an illustration, we consider an attractive 5-function potential in two dimensions, 
V{r) = — A2(5^(r). Although a trivial extension of the one- dimensional problem, it is rarely 
discussed^ in introductory courses, perhaps because the bound-state real-space wavefunc- 
tion is logarithmically divergent^ in the vicinity of the (5-function. The bound-state en- 
ergy Eb{X2, Uc) depends on the ultraviolet cutoff Uc and has a non-analytic dependence on 
the strength of the potential, Eij/Eq = — f/^ exp(— 47r£'o'2o/-^2)-— The momentum-space 
Schrodinger equation in this case is analytically tractable and provides a good test.— Because 
the Fourier transform of this potential is a constant, the Hamiltonian matrix becomes 

Hmn{,t} = ul6mn " ^"^^ TT^^ffl = — -f^nm(^)- (lO) 

The 5-function potential affects only i = sector of the Hilbert space because wavefunctions 
with i vanish at the position of the 5-function due to the centripetal barrier. We verify 
that there is a single bound state for an attractive potential (A2 > 0) and none for a repulsive 
potential (A2 < 0). Figure H] shows the magnitude of the bound-state energy -Efe(A2) as a 
function of 4:TrEQal/X2 for 2 < \2/EQal < 20. We use Am = 0.01 and two ultraviolet 
cutoffs Uc = {10,20}. At large values of \2/EQa^ ~ Uc-, the numerical results deviate from 
the expected straight-line behavior due to discretization. This deviation is systematically 
suppressed by reducing A-u. Note that for an attractive (5-potential in both one and two 
dimensions the bound-state wavefunction has the same functional form, tpbiv) oc l/(p^ + 
h'^K^). However, because the bound-state energy Eb2 oc exp(— I/A2) in two dimensions'^, in 
contrast to the bound-state energy Ehi oc in one dimension, ■i'^ the size of the wavefunction 



in momentum-space in two dimensions is much smaller than that in one dimension, hK2 = 
^j2m\E\)2 <^ = A/2m|i?fei|. We emphasize that the dependence of E^, on the cutoff Uc is 
a peculiar property of the weakly bound state in the two-dimensional 5-function potential 
and arises due to the absence of an energy scale in a problem characterized by m, A2). 
For a general potential, including the attractive Coulomb interaction Vir') = — e^/r, we 
numerically obtain multiple bound-states with energies that are independent of the cutoff.— 
For a central potential it is straightforward to extend this method to higher dimen- 
sions. In /^-dimensions we represent a vector using hyperspherical co-ordinates, p = 
{p,(j),6i,62, ... ,61)^2) where G [0,27r] and 6i G [0,7r].— The Hamiltonian is then block- 
diagonalized into blocks with different angular momenta. The effective potential in the 
block {£, iz) is obtained by performing an integral over angular variables, similar to that 
in Eq. iQ. It is not always possible to analytically carry out this integration. The re- 
sulting Hamiltonian matrix satisfies Hnm,{£,£z) = i'^m/unJ^^^Hmni^^^z), and the resulting 
eigenvectors ipa{uk) are orthogonal with respect to the D- dimensional inner product 

A ^ 

m^p) = = (« ^ p). (11) 

ZTT ^ — ^ 

fe=0 

The exploration of the Hamiltonian matrix Hmni^Az) in D > 2 dimensions emphasizes 
two important points. First, by explicit construction, it generates a set of matrices, each 
element of which appears non-Hermitian and still has a purely real eigenvalue spectrum. 
Second, it explicitly demonstrates that the notion of orthonormality and Hermiticity are inti- 
mately connected to the inner-product used to construct the Hilbert space of wavef unctions .- 

IV. CONCLUSIONS 

We have presented an approach to the real-space Schrodinger equation via the Hamil- 
tonian matrix in momentum-space that is obtained after a suitable discretization.- This 
method does not suffer from the instability associated with discretization of the real-space 
Schrodinger equation,-^ primarily because the kinetic energy term is diagonal in momentum 
space and, for most physical potentials, the amplitude Vppi for scattering from p to p' de- 
cays for large \p — p'\. Therefore, the elements of the matrix i/mn near the top-right and 
bottom-left corners are small. 

Our method is best suited for numerically investigating the energies and wavefunctions 
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of bound states that occur in a localized central potential V{r) with a finite Fourier trans- 
form V{q). Many well-known examples with confining potentials where all eigenstates are 
localized (an infinite quantum well or a simple harmonic oscillator) cannot be studied using 
our approach because the Fourier transform is ill-defined. However, as we have discussed 
in Sec. [ITl it is possible to explore the low-lying eigenstates of such a system by choosing 
parameters such that Vo/Ea = Vq/ {li? /2ma?) ^ 1. Such a deep well, as far as the low-lying 
eigenstates are concerned, can be treated as an infinite well. 

V. SUGGESTED PROBLEMS 

Problem 1. Obtain the bound-state spectra for the potential 

f-K,[l - (21x1/0)"] \x\<a/2 
V,{x) = (12) 

I = otherwise, 

where 77 > 0. Note that K,(x) represents a family of potentials that extrapolate from a linear 
(?7 = 1), a quadratic (?] = 2), to a quantum well {j] 00). Choose Vo/Ea = Vq/ {h? /2ma?) 3> 
1. Compare your results to the WKB approximation prediction Ek = —Ari[2k + ij^WCs+f?)^ 
where is a constant and the eigenenergies Ej^ are measured from the bottom of the 
potential well. 

Problem 2. Obtain the analog of Eq. ffTOj) in three dimensions and study the spectrum 
for a given cutoff Uc- Show that a bound state arises only when A3 > ~ and 
determine Asc. Contrast your results with those for a quantum well with depth Vq and size 
h/Pc = ao/U,. 
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Figure captions 



FIG. 1: (color online) (a) Dependence of the magnitude of the bound-state energy \E\,\ obtained 
from the matrix Hmn, Eq- ([5]), on the strength Ai of the one-dimensional attractive (5- function 
potential. The energy is in units of Eq = /2ma^ and Ai is in units of E^aQ. The numerical result 
(crosses) is in excellent agreement with the analytical result-'- \Ei)\/Eq = Af /4(£'o'^o)^ (dashed 
line), (b) Typical momentum-space wavefunctions for the bound state (top curve) and a positive- 
energy state (bottom curve) for Ai/E'oOo = 0.5. The momentum p is in units of ao/h and the 
wavefunction is in units of -^/oq. The width of the bound-state wavefunction "ipbip) is given by 
hn = mXi/h. The positive-energy wavefunction is, as expected, sharply localized in momentum 
space. 
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FIG. 2: (color online) Eigenenergies En of bound states in a deep quantum well Vo/Ea ^» 1 
obtained from Eq. ([5]). The bound-state energies En are positive because they are measured from 
the bottom of the well. The energies are in units of Ea = h'^/2ma'^ = Eo{ao/a)'^ instead of the 
customary unit Eq. The solid and the dashed lines represent results for Vo/Ea = (a/ao)'^ = 400 
and 625 respectively. The dotted line shows the analytical result for an infinite quantum well of 
width a, En = {n7r)'^Ea- 

FIG. 3: (color online) Magnitude of the ground-state energy \Eh\ for a Gaussian potential V{x) = 
— Vo exp(— x^/2a^) as a function of a for a fixed depth Vq/Eq = 1. The energies are in units of £"0, 
the length is in units of ao, and the momentum in the inset is in units of h/ao- The inset shows 
the ground-state momentum-space wavefunction in units of ^/ao■ As a increases, the ground-state 
wavefunction becomes increasingly localized in real space, and is reflected in the broadening of the 
momentum-space wavefunction. 

FIG. 4: (color online) Dependence of \Eh\ on the strength A2 of the attractive two-dimensional 
5-function potential for ultraviolet momentum cutoffs Uc = 10 (squares) and Uc = 20 (circles). 
The energy is in units of Eq and A2 is in units of Eoo^; An = 0.01 and 2 < A2/-Eo'^o — "^^^ 
solid lines represent the analytical result, ln{\Ef,\/ Eq) = —AirEQaQ/ X2 + ln([/^) for Uc = 10 (bottom 
curve) and Uc = 20 (top curve). The corresponding y-intercepts are In(lO^) ~ 4.6 and ln(20^) ~ 6. 
The discrepancy between the numerical and analytical results for large X2/EQaQ ~ Uc is decreased 
when Au is reduced. 
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